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CLASSICAL SCATTERING CALCULATIONS FOR 
DIATOMIC MOLECULES 
A General Procedure and Application 
to the Microwave Spectrum of 

Uri Mingelgrin 
ABSTRACT 

Many properties of gaseous systems such as 
electromagnetic absorption and emission, sound 
dispersion knd absorption, etc. may be elucidated 
if the nature of collisions between the particles in 
the system is understood. In this report, a proce- 
dure for the calculation of the classical trajectories 
of two interacting diatomic molecules will be 
described. The dynamics of the collision will be 
assumed to be that of two rigid rotors moving in a 
specified potential. The actual outcome of a 
representative sample of many trajectories at 298K 
was computed, and the use of these values at any 
temperature for calculations of various molecular 
properties will be described. Calculations performed 
for the microwave spectrum are given to 
demonstrate the use of the procedure described. 

Key words: Classical scattering, line shape, collision 
cross sections, O^microwave spectrum, O self 
broadening, spectrum foreign gas broadening. 

1. INTRODUCTION 

The elucidation of properties of an ensemble of particles in 
equilibrium requires the knowledge of ensemble averages, rather than 
the time development of the individual particles, over long periods of 
time. In addition, many properties of molecules in a gaseous system 
can be described by the impact approximation (Anderson, 1949) in the 
pressure range where the molecules are unperturbed except for the 
very short duration of binary collisions. To describe a trajectory of a 


molecule during such a binary collision, one needs to define the 
colliding molecule's initial conditions such as relative velocity, 
angular momenta, etc. , and then to solve the equations of motion 
assuming some intermolecular potential. One may, therefore, choose 
a statistical sample of initial conditions for the various possible colli- 
sions and then solve the proper equations of motion. Once a set of 
initial and final conditions for a sample of trajectories is given, cross 
sections and expectation values for various variables can be derived 
which can then be used for calculations of the various relaxation 
phenomena and other properties. For example, the microwave 
spectrum, the rotational Raman spectrum and sound absorption by 
were computed from the data derived from such a statistical set of 
collisions (to be published, also see section 3). 

Section 2 of this report will describe how a representative sample 
of collisions is calculated, while section 3 will demonstrate the use of 
such data in calculation of various properties. A description of a pro- 
gram which performed the trajectory calculations is given in section 
2. 3. 

2. THE CLASSICAL TRAJECTORY CALCULATIONS 

The collision between two diatomic molecules was assumed to be 
between two rigid rotators moving in a specified potential. This 
approximation will hold when: 
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a) the duration of practically all collisions is much longer than the 
vibrational period, so that (in analogy with the Born-Oppenheimer 
approximation) the vibrational motion can be approximately averaged in 
each intermolecular distance to give essentially two rigid rotors. 

b) the total energy associated with the rotational and relative motion is 
lower for most collisions than the energy difference between the 
vibrational levels. 


c) the energy differences between the rotational levels are significantly 
smaller than the energy difference between the vibrational levels. 

When conditions b and c hold, vibrational transitions may be neglected. 
All three conditions are met by O , N and many other diatomic 

Li Li 

molecules up to temperatures higher than room temperature. 

2. 1 Solution of the Equation of Motion 
To determine the trajectory of two colliding rotors it is necessary 
to solve the appropriate equations of motion. This can be done after 
the interaction potential between the two rotors is defined. For 
example, the potential chosen for the 0 -0 calculations consists of a 

Lu Ct 

4-point exponential repulsion and of an attractive, angle dependent 
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where V , a , C , T and A are potential parameters obtained 
o 

by fit to virial coefficient, spectroscopic data and theoretical calcula- 
tions (Mingelgrin, 1972). None of the parameters was adjusted to fit 
the microwave spectrum the calculation of which is described in section 
3. The values of the potential parameters used are: 

V = 48838 eV 
o 

a = 5. 3883 A -1 
C = 48.108 eVA 6 
r = . 229 
A = . 052 

The r.. are the distances between four centers of repulsion, two on 
ij 

each molecule. In this work the r. . were taken as the distances 

ij 

between the nuclei of one molecule and the nuclei of the other 

molecule. In addition, 0^ and 0^ are the angles between either 
molecule and the intermolecular axis, cc^-cp^ * s t ^ ie an g^ e between the 
projections of both molecules on a plane perpendicular to the inter- 
molecular axis and represents the second Legendre polynomial. 

The coordinate system chosen to describe the motion of the 
rotors is as follows. The relative motion of the two rotors is 
described in a Cartesian coordinate system with the initial vector of 
relative motion defining the Z axis. The axis perpendicular to it in 
the plane of motion is the Y axis. The X axis is thus perpendicular 
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to the plane of motion. This coordinate system was used and is 
described in detail by Karplus et al. (1965). The internal motion of the 
rotors is described in spherical polar coordinates. The axes defining 
the internal motion are parallel to those defining the relative motion. 
Using this coordinate system, Hamilton's equations of motion 
(Goldstein, 1950) are solved in the scattering program described below. 
The equations of motion are 


9H . 

dpj q j 


(2a) 


6H 

V' P ' 


(2b) 


where q^ are the specified coordinates and p^ are their conjugate 

momenta, q. and p. denote derivatives with respect to time. H is 
J J 

the Hamiltonian of the system. In terms of the chosen coordinate 
system, the Hamiltonian is defined as 

2 


3 p 2 , / . cp. 

H =£^ZiH p B. + -^7' + v(r) 

rrr r— ■ i \ i sin 9. 


(2c) 


3=1 


i=l 


where u, is the reduced mass for the relative motion; p. the conjugate 

momenta of the Cartesian relative coordinates; 9 and cp the angles 

l i ° 

defining the spherical polar coordinates of rotor i , cp. being the 
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azimuthal angle: I. is the moment of inertia of rotor i , and Y(r) 

the potential defined in equation 1. Finally, P , P are the 

0. cp. 

1 X 1 

appropriate conjugate momenta. 

The solutions of the above defined equations of motion were done 
by a variable step Nordsieck integration scheme (Nordsieck, 1962). 

The trajectory is considered complete when the intermolecular distance 
is larger than the initial intermolecular distance. The initial inter- 
molecular distance, in turn, is defined for each individual trajectory 
by the program scattering described below and is generally 20-26 a. u. 
The governing relation for determination of the initial intermolecular 
distance is: 

26 

R = p (.l) a ‘ U * (3) 

where R is the initial intermolecular distance and P is the magnitude 
of the initial relative momentum in atomic units. The change in the 
variables of the motion upon passing a fixed distance during the 
collision will have some inverse relation to the relative momentum. 

The actual relation (eq. 3) is an empirical one derived from test runs 
of a sample of 27 trajectories so that the difference between the value 

of any component of the linear and angular momenta at the chosen R 

_3 

and at 35 a. u. is less than 3 X 10 a.u. It is then assumed that the 
trajectory will approximate well a collision with R-^-« . 
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As is evident from equation 2c, there is a singularity in the 
Hamiltonian at 9^ = 0 . This singularity is not of a physical origin 
and thus can be removed for a given trajectory simply by changing the 
coordinate system. A 180° rotation of all the vectors defining the 
initial conditions about the line connecting the centers of mass of both 
rotors is performed when the calculations fail because of this apparent 
singularity. In this way an equivalent collision which will not pass 
through the previously encountered singularity is obtained. The 
internal Cartesian coordinates, the components of the angular momenta 
of the rotors and the Cartesian relative conjugate momenta are 
transformed in a straightforward way. The spherical coordinates and 
momenta used in the program are then calculated. The relative 
coordinates are not changed by the rotation. 

The program in its current version calculates, in addition to the 
variables of the motion, the classical rotational phase shifts for both 
rotors, defined as the angular change in the position of the rotor in the 
rotation plane during a collision (Gordon, 1966). This quantity is not 
well defined in non-instantaneous collisions. Before a collision the 
rotor rotates in some plane. After the collision is over, the rotor 
rotates again in some plane, different in general from the original one. 
In an instantaneous collision, the rotor will change its plane of rotation 
at a point where the initial and final rotational planes intersect. If no 
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phase change would have occurred in the collision, the rotor will rotate 
up to the instant of collision in the initial angular frequency and follow- 
ing that instant in the final angular frequency. At some time t after 
such a collision without a phase shift, the position of the rotor in the 
rotation plane is well defined. The angle at time t between the above 
defined position and the actual position of the rotor is the classical 
rotational phase shift. If the collision is not instantaneous, and during 
the collision both points of intersection of the initial and final planes of 
rotation are passed one or more times, the classical phase shift is not 
well defined. The procedure adapted for defining the classical phase 
shift is as follows. The time of closest approach in a collision is 
determined. Then, the two apparent phase shifts were calculated by 
assuming the rotor changed planes of rotation at both times of inter- 
section of the initial and final rotational planes nearest to the time of 
closest approach. These two apparent phase shifts are designated as 
and • The difference between the two points of intersection of 

the initial and rotational plane is, of course, tt radians. If we call the 
angles between the position of the rotor, if it had rotated up to the time 
of closest approach in the initial angular frequency and the points of 
intersection of the initial and final rotational planes, j3^ and (3^ > then 
the weighted average values of the various functions of the phase shift 
(F(v) ) are given by 
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F(v) = J + F(v 2 

2.2, Statistical Sampling and the Effect of Temperature 
To obtain a representative sample of trajectories, a procedure 
formulated by Conroy was used (Conroy, 1967); a Conroy case II for 
non-periodic variables with his parameters for 12 variables and 1861 
points was utilized. Conroy's procedure for multivariable sampling 
selects sampling points so that each variable is sampled uniformly 
along its range. For periodic variables one period is sampled. For 
non-periodic variables a transformation must be applied to let the range 
of the variable be 0 to 1 . For every variable a parameter k , which 
is the ratio between some odd integer smaller than the sample size and 
the number of points is defined. The sample size is chosen to be some 
large prime number. The values of a variable (x.) sampled may be 
selected by a number of schemes. In the case of "Conroy case II" the 
values selected are: 



x. = j • k - e . 
J x j 


j - 1,2, ...,•& 


(4) 


where & is the sample size and e. is an integer chosen so that x. 

J 3 

falls between -1, 1 . The absolute value of x. is then taken. Conroy 

J 

demonstrated that one set of points selected by such a procedure will 
optimize the sampling. The actual parameters k chosen were 
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selected by Conroy for various numbers of variables and sample sizes 

by trial and error to minimize the error in a test integral. The test 

integral was such as to maximize the error in integration. The 12 

variables chosen to completely define the initial conditions of a collision 

are: the absolute values of the angular momenta of both rotors; the 

impact parameter; the initial relative speed; the angles cp^ , cp^ , 9^, 

0 defined above; the conjugate momenta P and P ; and the signs 
2 9 1 6 2 

of P , P . The number of trajectories (1861) yields a sample of 
9 1 ^2 

3722 collision outcomes if both colliding rotors are identical, as in the 
case studied ( 0 -0 collisions). 

Li L 

As mentioned above, the initial values of the sampled variables 

had to be transformed so that the sampling range will be 0 to 1 . In 

the case of the angular momenta, a cutoff was selected so that the 

maximum rotational quantum number (J ) was 25; namely: 

max 

K = nV25(25+l)‘ , where K is the rotational angular momentum, 
max ° 

The initial rotational angular momenta allowed were only those which 
correspond to the odd rotational quantum numbers. The actual 
selection of the rotational angular momenta of the two rotors is as 
follows : 
let 

Q (J) = (2N+l)e _B N(N+1) /kT (5a) 

N=l, 3, . . . , J 
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and 


(J) Q(J ) 

max 


where B is the rotational constant. 

The range of P is 0 to 1 . The variable sampled by the 

( J) 

Conroy procedure is P . All points in the range 0-1 that fell 

w ) 

between P and P (or 0 and P(l) for J=1 ) were assigned to 

(J / 

the rotational angular momentum K = ikrrm . This procedure will 
sample J according to its Boltzmann distribution. Note that although 
the trajectories were classical, the selected initial angular momenta 
were properly quantized as only odd rotational quantum numbers are 
allowed for . 

For the sampling of the velocity and impact parameters at a fixed 
temperature, a procedure developed by R. G. Gordon (unpublished) was 


used. Consider a cross-section of the form 



dB < f(V, B) > G(B) G(V) 


( 6 ) 


where f(V, B) is a function of an individual trajectory. The function 
f(V, B) may be any function which depends on the trajectory. For 
example, table 1 gives a matrix whose elements are the corresponding 
functions f(V, B) needed for the calculation of the microwave 

spectrum at any density and temperature (see section 3). The variable 
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a in table 1 is the angle of reorientation of the rotational angular 
momentum in a collision and P is the probability of transfer in a 
particular collision from rotational state i to rotational state f . 
The angular brackets in equation 6 indicate average overall initial 
conditions other than relative velocity and the impact parameter. 
Here, V is the relative velocity, B the impact parameter. G(B) 
and G(V) are the probability functions for a collision to occur with 
the variables V and B for bimolecular collisions. 


G(B ) = 2ttB 


(7a) 


G(V) 




(7b) 


where V is the mean relative velocity, k is the Boltzmann constant, 
T the absolute temperature, and y the reduced mass for the collision. 

The transformation used to change the integration range from 
0 - <» to the desired range of 0 - 1 will now be described. Define the 
unnormalized weights 


-K (yV / 2kT) 

«B(V) = V e 


w(B) = (1+K1 B «A B M/ )" 1 (K2 b R^)" 1 


(8a) 


(8b) 
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where K^. , Kl^ , K2^ are dimensionless parameters and 

was selected to be the position of the minimum in the spherically 
averaged intermolecular potential. These parameters were selected so 
that < (b\ > - * s a PP rox i ma t e ly constant throughout the integration 
range, thus improving the sampling (see equation 11). The normalized 
weights are now defined as 



cu(X) 

U)(X) G(X) dX 


(9) 


where X represents V or B . Now if we use the transformation 

dA(X) = U(X) G(X) d(X) (10) 

the integration limits are 0 - 1 . We can now use Conroy's procedure 
to define values of A(V) and A(B) in the range 0-1 and then extract 
the appropriate values for V and B . Integral (6) is evaluated now 
from the sum 



< f(V,B) > 
U(V) U(B) 


(ID 


A sample calculation for deriving such an integral will be 
demonstrated in section 3. 
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The distribution functions for V and the rotational angular 
momenta are temperature dependent. Thus, the sampling described 
above is valid for one temperature only. However, over a wide range 
of temperatures it is possible to utilize the trajectory calculations for 
the reference temperature (298° K) by using appropriate conversions. 
Define 

W vb (T i ) = (U(B) U(V, T x ) f 1 (12) 


then 


W (T )U(V,T ) 

W (T ) = — — — 

VB V 2 ; U(V, T ) 

L 


(13) 


The value of a term in the sum (11) at given values of V and B , and 
temperature T is now 


W vb (T i) < f(V,B) > U(V, T x ) 
U(V,T 2 ) 


(14) 


In the Conroy procedure used, the values of V defining the terms of 
the sum of equation (11) were selected at fixed intervals in the function 
A(V) (equation (10) ) so that, at the reference temperature T^ , the 
integral (6) can be evaluated by the sum in equation (11) using the 
sample of V values directly. At an arbitrary temperature, however, 
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if we use the V values selected at T , we need to modify the sum 
for a proper estimate of the integral. For a large enough sample we 
have 

dA(V,T ) dA(V,T ) /dA(V,T A " 1 
AA(V ’ T 2 )s d\ T~ = — — 1 dV ) iA < V - T l> < 15 > 

From equation (10) we. obtain 

U(V,T )G(V,T ) 

iA < V ’ T 2» * PIV.T^V.T,) iA<V ’ T l» • < 16 > 

At the arbitrary temperature the necessary cross section 

can now be expressed as 

. V W VB (T 1 ) < f(V ’ B) > G(V,T 2 ) 

<T 2 ! " G(V,Tj) 

It is now necessary to properly^evaluate the average < f(V, B) > 
at an arbitrary temperature using the same sample of rotational 
angular momenta. 

This is a straightforward procedure. If we define the distribu- 
tion function of the angular momenta 


16 



G(J) = 


(2J t ne- B ° J(Jtl)/kT 

2 (2J + l)e- B ° J(J+1)/kT 

J=l, 3, . . . , Jmax 


(18) 


where J is the rotational quantum number of the perturbing molecule 
and B° the rotational constant, the average at an arbitrary tempera- 
ture will be obtained by evaluating the function f(V, B, J) at the J 

values selected for the distribution at the reference temperature 
and multiplying by the factor 


g(j,t 2 ) 

G(J, T x ) 


(19) 


Finally, at the arbitrary temperature the integral in equation 6 will be 
calculated from 



W VB (T 1 )G(V..T 2 )G(J.,T 2 )f(V i ,B.,J.) 

> LJ 

~ G(V.T 1 (G(J i ,T 1 ) 


( 20 ) 


where summation is over all trajectories with the same initial angular 
momentum for the absorbing molecule and N is the number of such 

*J 

trajectories. Note that the cross sections are averaged over the 
angular momenta of the perturbing molecule only, since in the case of 


the microwave spectrum of O and many other properties (see for 

Ui 
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example section 3), these are the necessary averages. Averages for 


both molecules' angular momenta can be treated similarly, simply by 
summing overall values of the rotational angular momenta of both 
molecules using the proper Boltzmann factors. For example, equation 
(20) will have to be multiplied by the factor 

G ( J 1 - T 2 ) 

G(J i r , T x ) 

x 

where J. indicate the rotational quantum number of the radiator and 
the sum will be over all trajectories. 

2. 3. Description of Program and Output^ 

The 4-body trajectory calculations described above are executed 
by a classical scattering program fashioned after a reactive scattering 
program written by K. Morokuma and L. Pedersen. The present 


A list of the program described below as well as the output for a 
complete set of trajectories at 298°K can be obtained on request from 
the ITS MM Wave Transmission Spectroscopy Section. (See Table 3 for 
sample output and text for output description. 
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program is written by the author for the CDC 6600 and 7600 computers. 
An average trajectory calculation takes 4 seconds on the CDC 6600. 
Reference will be made to the actual variables and subroutines' names 
to simplify the understanding of the program listing for those inter- 
ested in obtaining it. 

The following items will be described: 

a. A list of subroutines and their functions is given in table 2. 

b. The input to the program. A list of the input appears in table 4. 
The subprogram and statement numbers of the format statements 
according to which the input variables are read also appear in 
table 4. A sample input appears in figure 5. 

c. A description of the way the program executes a trajectory 
calculat ion. 

d. The program output. A sample of the punch output is given in 
table 3. A list of the printed output parameters is given in 
table 6. 

The program input consists of a set of 5 cards per set of 
trajectories followed by two cards per trajectory. A sample of the 
input is given in table 5. The input variables, definitions, and the 
location and statement number of the associated format statements are 
listed in table 4. The input includes some options regarding the form 
of output, the conditions for the set of trajectories and data which 
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TABLE 2 


The subroutines of the scattering program and their functions. 

Subroutine Entry Purpose 

Name 

ZINIT To define initial coordinates, momenta, 

and error and step size parameters for 
the Nordsieck integrator. 

" ZNITY To redefine initial conditions if a 

singularity is encountered in the 
trajectory path. 

ZDIFFE Nordsieck integrator for the solution of 

the equations of motion. 

ZEVAL Calculates time derivatives of coordi- 

nates and momenta at every step on 
trajectory path. 

ZQUNO To read input data for the set of trajec- 

tories to be calculated, as well as for 
each individual trajectory. To determine 
the initial angular momenta of the two 
rotors. 

ZOUPUT Defines all printed output. 

ZEVR Calculates rotational energy of rotors. 
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Continuation of Table 2 


ZVECPR 


XIT 


ZVIBR 


ZENERG 


Calculation initial and final components 

of the angular momenta vectors of 
rotors and relative motion, final 
velocity and kinetic energy. 

Special exit in case of program failure 

(e. g. due to defective input). 

Calculate initial magnitude of linear 

momentum of rotor. 

Defines value of potential energy at each 

integration step. 

ZNERGY Defines various constants and potential 
parameters. 
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TABLE 3: A SAMPLE OF TRAJECTORY OUTPUT AT REFERENCE TEMPERATURE 

2 9 8 ° K (SEE DEFINITION OF SYMBOLS IN TEXT SECTION 2.3). 



determine the initial conditions of individual trajectories. The condi- 
tions fixed for the set of trajectories are: The weights of the four 
atoms, the maximum integration time allowed per trajectory, a 
parameter defining the initial intermolecular distance, the number of 
trajectories per set, the temperature and the highest rotational 
quantum number allowed. The variables defining the initial conditions 
of a trajectory were listed in the text (section 2.2) and the actual input 
data per trajectory are given in input cards 6 and 7 as listed and 
defined in table 4. 

The basic features of the actual execution of a set of trajectory 
calculations will now be described. First, the necessary input for a set 
of trajectories is read. Next the various reduced masses (array W) are 
defined from the masses of the four atoms (array WP) in the main 
program. Subroutine ZINIT is then called. ZINIT calls in turn sub- 
routine ZQUNO where the array XU(N,II) is defined. This is the array 
of the functions P(J) defined in equation (5b) for both rotors. The 
determination of individual trajectories starts at this point. ZINIT is 
called again. It in turn calls ZQUNO where the input data for the 
individual trajectories are read and the quantum numbers J(array JAB) 
of the rotational angular momenta for both rotor s are determined. Upon 
returning to ZINIT the input data are used to define the initial coordi- 
nates and conjugate momenta (the arrays QI -- the array of the 
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The scattering program input 
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both rotors 
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DT defined as the maximum step size in the Nordsieck integration is presently computed in 
subroutine ZINIT for each trajectory. 


Sample Input 
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Table 6: The printed output 
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*This is the standard output for a successful trajectory. Error messages (e. g. for too long 
trajectories) will also be printed. 
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JAB (2) defined in line 15 for rotor 
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DIFF3 MAIN difference between the X, Y and Z 

DIFF4 components of the initial and final total 

DIFF5 angular momentum respectively x 10^ 

DIFF4 = the magnitude of the difference between the 
initial and final magnitude of the total 
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Cartesian coordinates of the system and PI-- the array of the spherical 
conjugate momenta for the internal motion and the Cartesian conjugate 
momenta for the relative motion). The parameters necessary for the 
Nordsieck integration (DT and ERROR) and the initial intermolecular 
distance ( V RS2 ) are also defined in ZINIT. As stated above, spherical 
polar coordinates are used for the internal motion of the rotors in 
solving the equations of motion. However, for some purposes, the 
Cartesian coordinates are preferable. This is the reason for defining 
the array QI above. The internal Cartesian coordinates will also be 
defined at each integration step. 

Finally, the various components of the energy and angular 
momenta are calculated after a call to subroutines ZVECPR. Upon 
returning to the main program the subroutine ZOUPUT is called to 
print the properties of the initial state. At this point the initial state of 
the collision is completely defined. The solution of the equations of 
motion follows. At every step of the integration the subroutine ZDIFFE 
is called. Here the Nordsieck scheme is executed. ZDIFFE calls 
subroutine ZEVAL, In ZEVAL, the time derivatives of the momenta 
and coordinates (the array DZ) are defined at every step. These 
derivatives are obtained through partial derivatives of the Hamiltonian 
with respect to the coordinates and momenta (see eq. (2) ). To get the 
above derivatives, the value of the potential at each step of the 
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integration is needed. The value of the potential is defined in subrou- 
tine ZENERG. ZENERG is called by ZEVAL every integration step. 

In ZENERG the potential (eq. (1) ) is defined in terms of the distances 
between the four atoms of the system (array R). Although the potential 
is in terms of the actual interatomic distances, the distance between the 
repulsion centers on a mol-ecule is not necessarily the bond length of 
the molecule but is rather an adjustable parameter. In other words, 
for calculation of the repulsive potential an apparent bond length is 
introduced. For the O calculations, however, the centers of 

w 

repulsion were assumed to be at the nuclear sites. 

ZENERG defines the potential of interaction and the user may 
replace ZENERG by another subroutine supplying any potential he 
chooses as long as it is in terms of the distances between atoms. 
ZENERG has to supply also the derivative of the potential with respect 
to the internuclear distances. The subroutine ZEVAL, before calling 
ZENERG, defines the internuclear distances from the present 
coordinates of the system. 

After ZDIFFE performs an integration step the main program 
checks if the trajectory has terminated. If not, the integration is 
continued. If the integration approaches the apparent singularity 
(eq. (2c) ), variable IDUD is set to 0 and the entry ZNITY in subroutine 
ZINIT is called. There, the rotation defined in section 2.1 is applied. 
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If the singularity is approached again the trajectory is cancelled. 

(About 1% of the trajectories had to be cancelled due to either two 
consecutive approaches to the singularity or any other reason, see 
below. ) 

After the integration of the equations of motion is completed, the 
output quantities are calculated and printed or punched. The final 
rotational angular momenta are defined in ZVECPR, as were the 
initial rotational angular momenta. The angles of reorientation of the 
angular momenta and the rotational phase shifts are defined in the 
main program. 

Every trajectory which was longer than 220000 a. u. (T.GT. 

220000) or had more than 1400 steps (IBKIN, GT. 1400) is back integrated 

to ascertain the success of the numerical forward integration. There 

is an option IDTL (defined in the input list table 4), which allows back 

integration and detailed print for every trajectory. In the production 

runs, back integration was performed only for trajectories which 

obeyed the above restrictions in trajectory time and number of steps. 

Back integration is performed by setting the time T=0 , reversing the 

signs of all final momenta and solving the equation of motion through 

calls to ZDIFFE. The back integration is set to be of the same length 

of time as the forward integration and it should then terminate at the 

initial coordinates of the forward integration. The final momenta should 
be the negative of the initial momenta of the forward integration. After the 
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final coordinates and momenta of the back integration are printed, the 
calculation of the next trajectory starts through a call to ZINIT for the 
definition of the initial conditions of the new trajectory and integration 
of the equations of motion as described for the previous trajectory. 

Tables 4 and 6 define the input and the printed output variables. 
Between them they define most key variables and thus no further list of 
variables is supplied. In addition, comments in the body of the program 
listing give further details. Some variables defined in the program are 
not used in the current version. These were defined for various 
purposes in the past and were not removed as they might be useful 
again in the future. Table 4 includes a few such variables marked with 
an asterisk. All calculations in the program are done in atomic units. 

The program output consists of printed and punched portions. 

The printed output gives details of the initial and final variables of the 
motion and the integration time. This printed output enables a check of 
the conservation of the components of the angular momenta as well as 
energy. In the printed output the coordinates printed are the Cartesian 
coordinates while the momenta are the spherical polar conjugate 
momenta for the internal motion and the Cartesian conjugate momenta 
for the relative motion. In addition, the results of the back integration 
are printed; A successful integration is determined by proper 
conservation of the constants of the motion and successful back 
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integration. A successful back integration was taken to mean that no 

coordinate or momentum, after back integration, differs by more than 
_2 

2 X 10 a. u. from its initial value. Any trajectory which failed in 
integration was repeated with tighter integration parameters. As 
mentioned above, about 1% of the trajectories were not accepted due to 
repeated failures of the integration from any cause. Table 6 describes 
the printed output. 

The punched output {a sample of which is reproduced in table 3) 

includes the initial and final properties of each collision necessary for 

the calculation of various properties and cross sections (e. g. Gordon, 

1968). The properties included in the punched output are: W (T ) -- 

VB 1 

defined in section 2. 2; cos a -- the cosine of the angular change in the 
orientation of the rotational angular momentum vector of a rotor in the 
collision; K. -- the initial angular momentum; -- the final angular 
momentum; cos v , sinv , cos 2v , sin2v , which are the average 
values of the sine and cosine of the classical rotational phase shift (v) 
and of 2v respectively, and V the relative initial velocity. The 
cosines and sines of both v and 2v are given since the values for 
these trigonometric functions are averages of two possible values as 
defined in the end of section 2.1. The averages for 2v cannot be 
derived directly from the averages for v (e. g. F(v) = cos 2v) 
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Since there are in the case of 0-0 collisions two identical 

L Li 

molecules per collision, each trajectory yields 2 sets of output. These 

appear in the two consecutive cards (or lines in table 3 ) with the same 

value for W. r _(T ) . In converting to different temperatures, it is the 
VB 1 

rotational angular momentum of the perturbing molecule in a set of two 
colliding molecules which corresponds to G(J.T) in equation (20). 

3. USE OF THE PROGRAM: O CALCULATIONS 

Cd 

The general usefulness of the data derived from the above pro- 
gram for a set of individual trajectories is in the fact that different 
properties can be calculated after the most complicated and computer - 
time consuming calculations, (namely the four-body trajectory calcula- 
tions) are done. Expressions for various cross sections necessary for 
the elucidation of a number of properties are given by Gordon (1968 and 
1967). The expressions in this section (equations (21)-(23)) are also 
derived in these last references. The calculations performed on the 
molecule microwave spectrum are described in the following to 
demonstrate one use of the data derived from the above program. 

The loss tangent can be defined as 

tan5( “ ) 'IkT Im £' - & ■ ive S ) " 1 - £■£= A-it < 21 > 

A is the absorption coefficient, K is the wavelength /2 tt , n is the 
number density of the absorbing molecule, p is the number density of 
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the perturbing molecule, k is the Boltzmann constant and uj is the 
angular frequency of the radiation times the unity matrix. T is the 
absolute temperature, d is a vector of all the magnetic dipole 
moment matrix elements, p is a diagonal matrix of the probabilities 
of the various molecular states, v is the mean relative velocity and 
cr is a relaxation matrix to be discussed below. All quantities except 
for ct are well known. The diagonal elements of cr have the following 
interpretation: 


-Impva. . = line shift 
11 


( 22 ) 


Repvcr.. = line width 
11 


for line i . For the microwave spectrum the line shift parameter 

equals zero as g is a real matrix. One can express a as 


cr = v _1 [ 2rrBd B < v(l - S) > (23) 

0 

where B is the impact parameter; v is the relative velocity; v is 
the mean relative velocity; S is a collisional transfer matrix for a 
single collision; and the brackets indicate average overall collisions 
with the same impact parameter. 
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A matrix element of the S matrix S. . defines the amplitude 

» ij 

transferred from spectral line j to spectral line i in a collision. 

(See Gordon, 1967, for the derivation and definition of the S matrix. 
Note, however, that the above S matrix does not define transitions 
between states as does the familiar scattering matrix. ) The microwave 
spectrum of in the frequency region discussed below is a multiplet 

spectrum, and thus every line arises from transitions between states, 
both with the same molecular rotational quantum number. Hence, each 
line can be assigned to a rotational level. The elements of the S 
matrix will be, therefore, a product of the probability of transition 
between the initial and final rotational level and another term defining 
the probability of transition between the various multiplet lines 
corresponding to each rotational level (Gordon, 1967). If we define 
f(V, B, J) in equation (20) to be 6_- S. . , equation (23) can be solved 

using equation (20) and the calculated set of trajectories (the trajecto- 
ries are completed as defined in section 2.1). A submatrix of 1 - S 
corresponding to the lines belonging to one rotational level of is 

reproduced in table 1 from Gordon (1967). The variable P^. in table 1 
is defined as probability of transfer in a particular collision from the 
initial rotational level i to some final rotational level f and a is 
defined as the angle of reorientation of the rotational angular momentum 
in a collision. 
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The final rotational angular momenta are not quantized due to the 
use of a classical scattering scheme. The scattering program 
described above does not quantize the final rotational angular momenta. 
For the sake of completeness the quantization procedure used for the 
calculations will now be described briefly. 

The final apparent quantum number J is defined through the 
relation 

K = £Vt(J+1) 


where is the actual final rotational angular momentum. Only odd 

rotational quantum numbers are allowed for O , making the difference 

Lt 

between the quantum numbers of neighboring levels 2 . After the 
collision, the molecule is assumed to have some probability of being in 
either of the two rotational levels with quantum numbers closest to J 
but not in any other rotational level. The probability of the molecule 
being after a collision in one of these two rotational levels with quantum 
number is taken as 



(24) 


where AJ is the difference between J and . Thus, equation (24) 

defines the variable P.. above. P„. is zero for all levels but the two 

fi fi 

with quantum numbers closest to J . The index i in P^. refers to 
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the initial rotational level in the collision which is well defined by the 
sampling procedure (see section 2. 2). 

A detailed discussion of the microwave computations and the 

comparison of experimental results to calculations are given by 
Mingelgrin (1972). Note that there are no adjustable parameters in this 
theory and that the potential parameters are derived from independent 

data (see section 2.1). 

2 

Some results are given in figures 1-4 and tables 7 and 8. Table 
7 gives the calculated line width parameters for the different lines and 
compares them to experimental results. There are two absorption 
lines, designated J+ and J- , per rotational level. The wide varia- 
tion in the experimental results makes the comparison with them diffi- 
cult. This points out the need for further study of the O spectrum at 

La 

low densities. Since the scattering calculations reported here are 
classical, the line width parameters for the lowest rotational level 
should have the largest error. 


A set of programs developed by R. G. Gordon's group at Harvard 
University designed to evaluate expressions of the form of equation 
(21) from results of scattering calculations, such as described above, 
was used to derive the following results. 
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TABLE 7 


The Line Widths of the Pure Microwave Spectrum 

(Temperature 298° - 300° K) 


Lines 

(J±) 

Half Width (MHz/mmHg) 



This _ . < 

Work B ’ ,S&G 

3 Z&M a 

A,S&G b 

A&G 6 

H&G C 

1 + 

2.3 8 


1.96 

2. 20 


1 - 

2.38 




1.97 

3 + 

2.36 

1.96 

1.71 

2.23 

2.07 

3 - 

2. 36 2. 09 

1.96 

1.92 



5 + 

2. 06 

1.56 


1.96 

1.80 

5 - 

2. 06 

1.60 

1.86 

1.99 


7 + 

2.11 

1.68 

2.05 

1.92 


7 - 

2.11 

1.70 


1.82 

2.0 1 

9 + 

1.96 

1.42 


1.93 


9 - 

1.96 

1,. 6 4 

1.97 

2. 00 

1.94 

11 + 

1.96 

1.60 

1.97 



1 1 - 

1.96 



1.97 


13 + 

1.92 0.87 

1.54 




1 3 - 

1.92 


1.88 

1.86 


1 5 + 

1.72 


1.77 



1 5 - 

1.72 



1.99 


17 + 

1.74 

1.50 




1 7 - 

1.74 


1.76 

1.82 


1 9 + 

1.59 


1.58 



19- 

1.59 


1.62 

1 . 91 


2 1 + 

1.61 





2 1 - 

1.61 0.83 


1.26 



23 + 

1.53 


1.26 



2 3 - 

1.53 


1.49 



a R. W. 

Zimmerer and M. 

Mizushima, 

Phys. Rev. 121, 

151 (1961). 

b R.S . 

Anderson, W. V. Smith, W. Gordy, Phys. 

Rev. 87, 

561 (1952). 

c R. M. 

Hill and W. Gordy, 

Phys. Rev. 

93, 1019 (1954). 



d J. H. Burkhalter, R. S. Anderson, W. V. Smith, W. Gordy, Phys. Rev. 79 , 


651 (1950). 

e J. O. Artman, Absorption of Microwaves by Oxygen in the 

Millimeter Wavelength Region (Columbia Radiation Laboratory, 1963); 
J. O. Artman and J. P. Gordon, Phys. Rev. 96 , 1237 (1954). 
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Selected Values for the Pure O Absorption Coefficient at 298°K about the O 9+ line (db/km) 
















































































In figure 1, some experimental and calculated spectra are 
compared. Maryott and Birnbaum (I960) pointed out the problems 
associated with fitting a theoretical line shape generated by addition of 
Lorentzian lines with empirical parameters to higher pressure data 
(approximately 10 atmospheres and up). Despite the fact that no 
parameters were adjusted to fit the data, figure 1 demonstrates the good 
fit achieved by the present computations. This fit of calculation to 
experiment improves with decreasing pressure. 

Figure 2 compares calculations to experimental determinations 
of the absorption in an O^-Ar mixture. Here we have to 

calculate two kinds of interactions, 0 -0 and O -Ar . The O -Ar 

Cd Ld Ld Ld 

collisions were determined by an available rigid rotor-sphere scatter- 
ing program. The O -Ar calculations are described by Mingelgrin 

Cd 

(1972). The 0 -0 were determined with the help of the present 

Ca Cd 

program. The O^-Ar calculations, although not completely derived 
by the program described above, were introduced to demonstrate the 
possibility of extending the calculations to gas mixtures such as the 
atmosphere. O^-N^ an< ^ °^ er O^-diatomic calculations can be made 
with the present program after very minor modifications. Then, cross 
sections for the O^ microwave spectrum in gas mixtures can be 
derived. 


46 


The scattering calculations (equation (21)), using a complete 
relaxation matrix, introduce "interference" between lines. Namely, 
at densities where the individual lines overlap, the total line is 
narrower than predicted from Lorentzian lines addition. This 
"interference" should increase with pressure. Thus, the scattering 
calculations should give higher attenuation at the center of the band 
(~60 GHz ) and lower attenuation at the wings. At low densities, where 
the individual lines are resolved, the two procedures should be equiva- 
lent. The equivalence of the two procedures at the low densities where 
the lines are resolved arises from the fact that the off diagonal elements 
in the a matrix approach zero and equation (21) becomes equivalent to 
addition of Lorentzian lines. 

Figures 3 and 4 compare results of the present calculations for 
the pure spectrum to results of addition of Lorentzian lines with 

the same line width parameters. The linear addition spectrum was 
computed by W. M. Welch, I. T. S. The densities are all in the range 
of densities existing in the atmosphere. It is evident that there is a 
significant difference between the two calculations. This difference is 
gradually declining as the pressure declines and the spectrum 
approaches the resolved spectrum. In figure 4 a comparison 

between the present calculations and addition of Lorentzian lines at 
.13 Atmospheres (100 torr) and 298° K is given. Even at this density the 
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CALCULATED 



48 


I IS THE ACTUAL RANGE OF EXPERIMENTAL DETERMINATIONS, 
DETERMINATIONS PER FREQUENCY -PRES SURE POINT 



0.80 
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EXP. VS. CALCULATED 0 2 SPECTRUM IN 0 2 -Ar MIXTURES (298°K) 

I IS THE ACTUAL RANGE OF EXPERIMENTAL DETERMINATIONS, 

7 DETERMINATIONS PER FREQUENCY -PRESSURE POINTS. 



two calculations differ. However, the difference is much smaller, as 


expected, than the difference at the higher pressures (see figure 3). 
Figure 4 demonstrates another fact of importance for applications. At 
densities where the spikes of the different lines are still prominent, the 
two calculations do not differ by a similar amount at neighboring 
frequencies. In general, the difference is larger at the peaks of the 
lines and the troughs between the lines and smaller in the intermediate 
frequencies. 

The fit to experiment demonstrated in figure 1, especially at the 
lower density (3. 74 Atmospheres), and the difference between the two 
methods of calculation described in figure 3 demonstrate the need for 
this more complex procedure for the purpose of determination of the 
spectrum in the atmosphere in accuracies sufficient for practical 
applications. One has to await experimental results for pure 
spectra in the pressure range of interest to unequivocally state that the 
scattering calculations represented in figure 3 give the good fit 
presumed. (Such measurements are presently being undertaken by the 
millimeter wave transmission spectroscopy section of ITS. ) However, 
the fit in figure 1 and its improvement with decreased pressure suggest 
that the procedure does yield accurate spectra in the atmospheric 
pressure range. 
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COMPARISON OF COMPUTATIONS -SCATTERING CALCULATIONS 
VS. ADDITION OF LORENTZI AN LINES. (PURE 0 9 AT 298°K) 





It may seem that the equivalent of the scattering procedure can 
be achieved by adding Lorentzians with pressure dependent shift and 
width parameters. This is exactly so, except that the procedure 
described here does not rely on an empirical fit. As such, the method 
is universal and enables calculations of complete frequency or density 
profiles. Empirical fit of pressure dependent width and shift para- 
meters will depend on the accuracy of available measurements and will 
be applicable only to a small range of conditions as determined by the 
experimental points to which the parameters are fitted. 

Finally, the dispersion spectrum of at the microwave region 
can be easily computed using the same information and programs used 
for computation of the absorption spectrum. The r-eal part of the 
refractive index at frequency uj (n(o))) can be obtained from the relation 

n(to) -n(o) = ■ Red(w-^ Q -ivpo) 1 Pd (25) 

where all symbols have the same meaning as in equation (21). Equation 
(25) is discussed more fully in a forthcoming paper by the author. This 
equation may not apply at sufficiently high density, but it certainly 
holds for all pressures around and below 1 Atmosphere at 298°K. Some 
results of dispersion calculations are given in figures 5 and 6. 
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298 0 K . (n(u) - n ( 0) , WHERE n(u>) IS THE^REAL PART OF THE REFRACTIVE 

INDEX AT FREQUENCY w. ) 





CD 


gOl X [ (0 ) U - (m)U] 
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-80 -60 -40 -20 0 20 40 60 80 

FREQUENCY (MHz FROM CENTER OF 9+ LINE) 

THE CALCULATED DISPERSION SPECTRUM OF 0~ ABOUT THE 9+ LINE AT VARIOUS 
DENSITIES AT 298°K. (n(w) - n(0), WHERE^n(w) IS THE REAL PART OF THE 
REFRACTIVE INDEX AT FREQUENCY a>.) 



Zeeman splitting was neglected in all the computations. This, 
however, will affect the calculations at 298°K and pressures above 5 torr 
only negligibly for magnetic field strengths below or around . 5 Gauss 
which is approximately the strength of the earth's magnetic field 
(Liebe and Welch, 1972). 
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